## Make the plot of raw pairwise differences in social cohesion

library(here)
library(tidyverse)
library(viridis) ## for colorblind color palettes
library(sp)
library(robustbase)
library(ggExtra)

load(here("Design", "matches_anyDA_new.rda"))
load(here("Data", "wrkdatOwnMap_anyDA_new.rda"))
source(here("Design", "nonbimatchingfunctions.R"))

wdat0 <- inner_join(wrkdatOwnMap_new, matches_anyDA_new)

## Add information about which unit is ranked higher versus lower within pair
wdat0 <- wdat0 %>%
  group_by(pair) %>%
  mutate(
    perc_more = rank(vm.community.subj) - 1,
    social_capital01_rank = rank(social.capital) - 1,
    community_resp01_rank = rank(community.resp01) - 1
  ) %>%
  ungroup()

## Make data at the level of the pair to illustrate the pairwise differences
paired_data <- wdat0 %>%
  filter(perc_more != .5 & !is.na(pair)) %>%
  droplevels() %>%
  group_by(pair) %>%
  summarize(
    social_cohesion_diff = social.capital01[perc_more == 1] - social.capital01[perc_more == 0],
    perceptions_diff = vm.community.norm2[perc_more == 1] - vm.community.norm2[perc_more == 0]
  ) %>%
  ungroup() %>%
  na.omit()

colors <- viridis(3)
names(colors) <- c("OLS", "Outlier Robust Linear", "Outlier Robust Loess")

g_pairs_sc <- ggplot(data = paired_data, aes(x = perceptions_diff, y = social_cohesion_diff)) +
  geom_point(alpha = .5, color = "gray") +
  geom_hline(yintercept = 0) +
  geom_smooth(
    method = "loess", se = FALSE,
    method.args = list(family = "symmetric"),
    aes(col = "Outlier Robust Loess")
  ) +
  geom_smooth(se = FALSE, method = "lm", aes(color = "OLS")) +
  geom_smooth(
    se = FALSE, method = "lmrob", aes(color = "Outlier Robust Linear"),
    method.args = list(setting = "KS2014")
  ) +
  xlab("Difference in Perceptions (high vs low)") +
  ylab("Difference in Social Cohesion (high perceiver versus low perceiver)") +
  scale_color_manual(name = "", values = colors) +
  theme_classic() +
  theme(legend.position = "inside", legend.position.inside = c(.8, .8))
g_pairs_sc

g_pairs_w_dens <- ggMarginal(g_pairs_sc, type = "densigram", margins = "both")

g_pairs_w_dens

ggsave(g_pairs_w_dens, file = here("Figures_Tables", "pairwise_plot_social_cohesion_anyDA_new.pdf"), width = 5, height = 5)
